#######################################################################################################
## Replication Script                                                                                ##    
## State Capacity, Refugee Reception, and Attitudes towards Refugees                                 ## 
## Study 3: Further Evidence of the Positive Effects of Expanding Reception Capacity                 ##  
#######################################################################################################

#######################################################################################################
#### packages ####
library(rio)
library(ggplot2)
library(scales)
library(patchwork)
library(modelsummary)
library(marginaleffects)
library(estimatr)
library(mediation)

set.seed(0987)

#######################################################################################################

#######################################################################################################
#### sessionInfo() ####

# R version 4.4.0 (2024-04-24)
# Platform: aarch64-apple-darwin20
# Running under: macOS Sonoma 14.5

# Matrix products: default
# BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib  
# LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

#locale:
#  [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

# time zone: America/Chicago
# tzcode source: internal

# attached base packages:
#  [1] stats     graphics  grDevices utils     datasets  methods   base   

# other attached packages:
# [1] mediation_4.5.1        sandwich_3.1-1         mvtnorm_1.2-5          Matrix_1.7-0           MASS_7.3-61           
# [6] estimatr_1.0.4         marginaleffects_0.21.0 modelsummary_2.1.1     patchwork_1.2.0        scales_1.4.0          
# [11] ggplot2_3.5.2          rio_1.2.0          

# loaded via a namespace (and not attached):
# [1] tidyselect_1.2.1    dplyr_1.1.4         farver_2.1.2        R.utils_2.12.3      fastmap_1.2.0       bayestestR_0.15.2  
# [7] digest_0.6.36       rpart_4.1.23        estimability_1.5.1  lifecycle_1.0.4     cluster_2.1.6       magrittr_2.0.3     
# [13] compiler_4.4.0      rlang_1.1.6         Hmisc_5.1-3         tools_4.4.0         data.table_1.16.0   knitr_1.48        
# [19] labeling_0.4.3      texreg_1.39.4       htmlwidgets_1.6.4   RColorBrewer_1.1-3  abind_1.4-5         tinytable_0.3.0  
# [25] purrr_1.0.2         withr_3.0.2         foreign_0.8-87      R.oo_1.26.0         datawizard_1.0.0    nnet_7.3-19      
# [31] grid_4.4.0          fansi_1.0.6         xtable_1.8-4        colorspace_2.1-0    future_1.33.2       emmeans_1.10.3   
# [37] globals_0.16.3      insight_1.0.2       cli_3.6.5           rmarkdown_2.27      generics_0.1.3      performance_0.12.2
# [43] rstudioapi_0.16.0   future.apply_1.11.2 httr_1.4.7          parameters_0.24.1   minqa_1.2.7         stringr_1.5.1    
# [49] splines_4.4.0       parallel_4.4.0      effectsize_0.8.9    base64enc_0.1-3     vctrs_0.6.5         boot_1.3-30      
# [55] carData_3.0-5       car_3.1-2           Formula_1.2-5       htmlTable_2.4.3     listenv_0.9.1       tidyr_1.3.1      
# [61] glue_1.8.0          parallelly_1.37.1   nloptr_2.1.1        codetools_0.2-20    stringi_1.8.4       gtable_0.3.6     
# [67] tables_0.9.28       lme4_1.1-35.5       tibble_3.3.0        pillar_1.11.0       htmltools_0.5.8.1   R6_2.6.1          
# [73] evaluate_0.24.0     lpSolve_5.6.23      lattice_0.22-6      R.methodsS3_1.8.2   backports_1.5.0     broom_1.0.6      
# [79] Rcpp_1.0.13         coda_0.19-4.1       gridExtra_2.3       nlme_3.1-165        checkmate_2.3.1     xfun_0.46        
# [85] zoo_1.8-14          pkgconfig_2.0.3   
#######################################################################################################

#######################################################################################################
#### Study 3: Table A13: Sample Descriptives and Randomization ####

data <- rio::import("Hlatky_UKR_Ref_Study3_Exp_Data.csv")

datasummary_balance(sex+age+edu+reg+size+income+lib_con~treated, data=data, output="latex", starts=c('*' = .05), fmt = fmt_decimal(digits = 3, pdigits = 3))

#######################################################################################################

#######################################################################################################
#### Study 3: Figure 9: Average Treatment Effects of EU Financial Support ####

data$ukr_therm <- data$ukr_therm * 0.1

modela <- lm_robust(ukr_therm~treated, data=data)
summary(modela)

modelb <- lm_robust(burden_benefit~treated, data=data)
summary(modelb)

modelc <- lm_robust(ukr_accept~treated, data=data)
summary(modelc)

modeld <- lm_robust(eu_feel~treated, data=data)
summary(modeld)

exp_models <- list("EU Attitudes"=modeld, "UKR. Refugees Feeling Thermometer" =modela, "UKR. Refugees Burden/Benefit"=modelb, "Accept Fewer/More UKR. Refugees"=modelc)

exp_model_plot <- modelplot(exp_models, size=.1,  coef_omit = c(-2), conf_level=.95, color="black", facet = TRUE, coef_rename = "", ) + geom_vline(xintercept = 0, linetype = 'dashed') + theme_minimal() + xlab("Average Treatment Effect") +  theme(axis.text.x = element_text(size = 8), axis.text.y = element_text(size = 8), axis.title.x = element_text(size = 8)) + xlim(-0.5,1)

exp_model_plot

# ggsave("Fig9_Exp_Results.jpeg", plot = exp_model_plot, width=6 ,height = 4, dpi=600)

#######################################################################################################
#### Table A14: ATE Estimation, Figure 9 Main Text ####

exp_models_short_label <- list("(EU)"=modeld, "(UKR. Ref. Therm.)" =modela, "(Burden/Benefit)"=modelb, "(Accept UKR Ref.)"=modelc)


cm_table_models <- c("treated"="Treatment")

modelsummary(exp_models_short_label, output = "latex", coef_map = cm_table_models, gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05))

#######################################################################################################

#######################################################################################################
#### Study 3: Figure 10: Conditional Average Treatment Effects of EU Financial Support across Individual Ideology ####

data$lib_con_binary <- ifelse(data$lib_con<=5,0,1)

data$lib_con_binary <- factor(data$lib_con_binary, 
                       levels = c(0, 1), 
                       labels = c("Liberal", "Conservative"))

model1 <- lm_robust(ukr_therm~treated*as.factor(lib_con_binary), data=data)
summary(model1)

lib_con_therm_cate_plot <- plot_slopes(
  model = model1,
  variable = "treated",
  by = "lib_con_binary",  
  conf_level = 0.95,
  gray = TRUE
) +
  scale_x_discrete(
    labels = c("", "")
  ) +
  scale_y_continuous(
    limits = c(-1.2, 1.2),
    breaks = seq(-1, 1, by = 0.5),
    labels = label_number(accuracy = 0.1)) +
  theme_minimal() +
  xlab("") +  
  ylab("Marginal Effect of Treatment") +
  labs(title = "UKR. Refugees Feeling Thermometer") +
  theme(
    plot.title = element_text(
      size = 9,         
      hjust = 0.5,      
      face = "plain"    
    ),
    axis.text.x = element_text(size = 7),  
    axis.text.y = element_text(size = 6.5),
    axis.title.x = element_blank(),  
    axis.title.y = element_text(size = 8),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    legend.position = "bottom",
    legend.title = element_blank(),
    plot.margin = margin(5, 5, 0, 5)
  ) + geom_hline(yintercept = 0, linetype = "dashed")

hist_plot_grouped <- ggplot(data, aes(x = lib_con_binary, fill = as.factor(treated))) +
  geom_bar(alpha = 0.7,width = 0.2, position = position_dodge(width = 0.2)) + 
  scale_x_discrete(labels = c("Liberals", "Conservatives")) +
  scale_y_continuous(limits = c(0, 425), breaks = seq(0, 425, by = 200)) +
  scale_fill_manual(
    values = c("1" = "black", "2" = "gray40"),
    labels = c("1" = "Control", "2" = "Treatment"),
    name = ""
  ) +
  theme_minimal() +
  xlab("") +
  ylab("Count") +
  theme(
    axis.text.x = element_text(size = 7),
    axis.text.y = element_text(size = 5),
    axis.title.x = element_text(size = 7),
    axis.title.y = element_text(size = 7),
    panel.grid.minor.y = element_blank(), 
    panel.grid.minor.x = element_blank(),
    legend.position = "bottom",
    legend.text = element_text(size = 6),
    legend.key.size = unit(0.3, "cm"),
    plot.margin = margin(0, 5, 5, 5)
  )

model2 <- lm_robust(burden_benefit~treated*as.factor(lib_con_binary), data=data)
summary(model2)



lib_con_burden_cate_plot <- plot_slopes(
  model = model2,
  variable = "treated",
  by = "lib_con_binary",  
  conf_level = 0.95,
  gray = TRUE
) +
  scale_x_discrete(
    labels = c("", "")
  ) +
  scale_y_continuous(
    limits = c(-1.2, 1.2),
    breaks = seq(-1, 1, by = 0.5),
    labels = label_number(accuracy = 0.1)) +
  theme_minimal() +
  xlab("") +  
  ylab("") +
  labs(title = "UKR. Refugees Burden/Benefit") +
  theme(
    plot.title = element_text(
      size = 9,         
      hjust = 0.5,      
      face = "plain"    
    ),
    axis.text.x = element_text(size = 7), 
    axis.text.y = element_text(size = 6.5),
    axis.title.x = element_blank(),  
    axis.title.y = element_text(size = 8),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    legend.position = "bottom",
    legend.title = element_blank(),
    plot.margin = margin(5, 5, 0, 5)
  ) + geom_hline(yintercept = 0, linetype = "dashed")






model3 <- lm_robust(ukr_accept~treated*as.factor(lib_con_binary), data=data)
summary(model3)

lib_con_accept_cate_plot <- plot_slopes(
  model = model3,
  variable = "treated",
  by = "lib_con_binary",  
  conf_level = 0.95,
  gray = TRUE
) +
  scale_x_discrete(
    labels = c("", "")
  ) +
  scale_y_continuous(
    limits = c(-1.2, 1.2),
    breaks = seq(-1, 1, by = 0.5),
    labels = label_number(accuracy = 0.1)) +
  theme_minimal() +
  xlab("") +  
  ylab("Marginal Effect of Treatment") +
  labs(title = "Accept Fewer/More Ukr. Refugees") +
  theme(
    plot.title = element_text(
      size = 9,          
      hjust = 0.5,       
      face = "plain"     
    ),
    axis.text.x = element_text(size = 7),  
    axis.text.y = element_text(size = 6.5),
    axis.title.x = element_blank(),  
    axis.title.y = element_text(size = 8),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    legend.position = "bottom",
    legend.title = element_blank(),
    plot.margin = margin(5, 5, 0, 5)
  ) + geom_hline(yintercept = 0, linetype = "dashed")


model4 <- lm_robust(eu_feel~treated*as.factor(lib_con_binary), data=data)
summary(model4)


lib_con_eu_cate_plot <- plot_slopes(
  model = model4,
  variable = "treated",
  by = "lib_con_binary",  
  conf_level = 0.95,
  gray = TRUE
) +
  scale_x_discrete(
    labels = c("", "")
  ) +
  scale_y_continuous(
    limits = c(-1.2, 1.2),
    breaks = seq(-1, 1, by = 0.5),
    labels = label_number(accuracy = 0.1)) +
  theme_minimal() +
  xlab("") +  
  ylab("") +
  labs(title = "EU Attitudes") +
  theme(
    plot.title = element_text(
      size = 9,          # make it small
      hjust = 0.5,       # center it
      face = "plain"     # plain text (no bold)
    ),
    axis.text.x = element_text(size = 7),  # now show your two labels
    axis.text.y = element_text(size = 6.5),
    axis.title.x = element_blank(),  
    axis.title.y = element_text(size = 8),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    legend.position = "bottom",
    legend.title = element_blank(),
    plot.margin = margin(5, 5, 0, 5)
  ) + geom_hline(yintercept = 0, linetype = "dashed")

top_plots_grid <- wrap_plots(
  lib_con_accept_cate_plot,
  lib_con_burden_cate_plot,
  lib_con_therm_cate_plot,
  lib_con_eu_cate_plot,
  ncol = 2,
  nrow = 2
)

bottom_plots <- hist_plot_grouped | (hist_plot_grouped + theme(legend.position = "none"))


combined_cate_plot <- top_plots_grid / bottom_plots +
  plot_layout(heights = c(6, 0.5), guides = "collect") & 
  theme(legend.position = "bottom")

combined_cate_plot

# ggsave("Fig10_LibCon_Cate.jpeg", plot = combined_cate_plot, width=6 ,height = 5, dpi=600)

#######################################################################################################


#######################################################################################################
#### Study 3: Table A15: CATE Estimation Figure 10 Main Text ####

exp_cate_models_short_label <- list("(EU)"=model4, "(UKR. Ref. Therm.)" =model1, "(Burden/Benefit)"=model2, "(Accept UKR Ref.)"=model3)

cate_coef_map <- c("treated"="Treatment", "as.factor(lib_con_binary)Conservative"="Conservative", "treated:as.factor(lib_con_binary)Conservative"="Treatment*Conservative")

modelsummary(exp_cate_models_short_label, output = "latex", coef_map = cate_coef_map, gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05))

#######################################################################################################
#### Study 3: Figure A7: Conditional Average Treatment Effects of EU Financial Support across Individual Ideology (Continuous)  ####

model1_c <- lm_robust(ukr_therm~treated*lib_con, data=data)
summary(model1_c)

lib_con_therm_cate_plot <- plot_slopes(
  model = model1_c,
  variable = "treated",
  by = "lib_con",  
  conf_level = 0.95,
  gray = TRUE
) +
  scale_x_discrete(
    labels = c("", "")
  ) +
  scale_y_continuous(
    limits = c(-2, 2),
    breaks = seq(-2, 2, by = 0.5),
    labels = label_number(accuracy = 0.1)) +
  theme_minimal() +
  xlab("") +  
  ylab("Marginal Effect of Treatment") +
  labs(title = "UKR. Refugees Feeling Thermometer") +
  theme(
    plot.title = element_text(
      size = 9,         
      hjust = 0.5,      
      face = "plain"    
    ),
    axis.text.x = element_text(size = 7),  
    axis.text.y = element_text(size = 6.5),
    axis.title.x = element_blank(),  
    axis.title.y = element_text(size = 8),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    legend.position = "bottom",
    legend.title = element_blank(),
    plot.margin = margin(5, 5, 0, 5)
  ) + geom_hline(yintercept = 0, linetype = "dashed")

hist_plot_grouped <- ggplot(data, aes(x = lib_con, fill = as.factor(treated))) +
  geom_bar(alpha = 0.7,width = 0.2, position = position_dodge(width = 0.2)) + 
  scale_x_continuous(limits = c(-0.5, 10.5),breaks = seq(0, 10, by = 1)) +
  scale_y_continuous(limits = c(0, 200), breaks = seq(0, 200, by = 50)) +
  scale_fill_manual(
    values = c("1" = "black", "2" = "gray40"),
    labels = c("1" = "Control", "2" = "Treatment"),
    name = ""
  ) +
  theme_minimal() +
  xlab("") +
  ylab("Count") +
  theme(
    axis.text.x = element_text(size = 7),
    axis.text.y = element_text(size = 5),
    axis.title.x = element_text(size = 7),
    axis.title.y = element_text(size = 7),
    panel.grid.minor.y = element_blank(), 
    panel.grid.minor.x = element_blank(),
    legend.position = "bottom",
    legend.text = element_text(size = 6),
    legend.key.size = unit(0.3, "cm"),
    plot.margin = margin(0, 5, 5, 5)
  )

model2_c <- lm_robust(burden_benefit~treated*lib_con, data=data)
summary(model2_c)



lib_con_burden_cate_plot <- plot_slopes(
  model = model2_c,
  variable = "treated",
  by = "lib_con",  
  conf_level = 0.95,
  gray = TRUE
) +
  scale_x_discrete(
    labels = c("", "")
  ) +
  scale_y_continuous(
    limits = c(-2, 2),
    breaks = seq(-2, 2, by = 0.5),
    labels = label_number(accuracy = 0.1)) +
  theme_minimal() +
  xlab("") +  
  ylab("") +
  labs(title = "UKR. Refugees Burden/Benefit") +
  theme(
    plot.title = element_text(
      size = 9,         
      hjust = 0.5,      
      face = "plain"    
    ),
    axis.text.x = element_text(size = 7), 
    axis.text.y = element_text(size = 6.5),
    axis.title.x = element_blank(),  
    axis.title.y = element_text(size = 8),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    legend.position = "bottom",
    legend.title = element_blank(),
    plot.margin = margin(5, 5, 0, 5)
  ) + geom_hline(yintercept = 0, linetype = "dashed")






model3_c <- lm_robust(ukr_accept~treated*lib_con, data=data)
summary(model3_c)

lib_con_accept_cate_plot <- plot_slopes(
  model = model3_c,
  variable = "treated",
  by = "lib_con",  
  conf_level = 0.95,
  gray = TRUE
) +
  scale_x_discrete(
    labels = c("", "")
  ) +
  scale_y_continuous(
    limits = c(-2, 2),
    breaks = seq(-2, 2, by = 0.5),
    labels = label_number(accuracy = 0.1)) +
  theme_minimal() +
  xlab("") +  
  ylab("Marginal Effect of Treatment") +
  labs(title = "Accept Fewer/More Ukr. Refugees") +
  theme(
    plot.title = element_text(
      size = 9,          
      hjust = 0.5,       
      face = "plain"     
    ),
    axis.text.x = element_text(size = 7),  
    axis.text.y = element_text(size = 6.5),
    axis.title.x = element_blank(),  
    axis.title.y = element_text(size = 8),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    legend.position = "bottom",
    legend.title = element_blank(),
    plot.margin = margin(5, 5, 0, 5)
  ) + geom_hline(yintercept = 0, linetype = "dashed")


model4_c <- lm_robust(eu_feel~treated*lib_con, data=data)
summary(model4_c)


lib_con_eu_cate_plot <- plot_slopes(
  model = model4_c,
  variable = "treated",
  by = "lib_con",  
  conf_level = 0.95,
  gray = TRUE
) +
  scale_x_discrete(
    labels = c("", "")
  ) +
  scale_y_continuous(
    limits = c(-2, 2),
    breaks = seq(-2, 2, by = 0.5),
    labels = label_number(accuracy = 0.1)) +
  theme_minimal() +
  xlab("") +  
  ylab("") +
  labs(title = "EU Attitudes") +
  theme(
    plot.title = element_text(
      size = 9,          # make it small
      hjust = 0.5,       # center it
      face = "plain"     # plain text (no bold)
    ),
    axis.text.x = element_text(size = 7),  # now show your two labels
    axis.text.y = element_text(size = 6.5),
    axis.title.x = element_blank(),  
    axis.title.y = element_text(size = 8),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    legend.position = "bottom",
    legend.title = element_blank(),
    plot.margin = margin(5, 5, 0, 5)
  ) + geom_hline(yintercept = 0, linetype = "dashed")




top_plots_grid <- wrap_plots(
  lib_con_accept_cate_plot,
  lib_con_burden_cate_plot,
  lib_con_therm_cate_plot,
  lib_con_eu_cate_plot,
  ncol = 2,
  nrow = 2
)

bottom_plots <- hist_plot_grouped | (hist_plot_grouped + theme(legend.position = "none"))


combined_cate_plot <- top_plots_grid / bottom_plots +
  plot_layout(heights = c(6, 0.5), guides = "collect") & 
  theme(legend.position = "bottom")

combined_cate_plot

# ggsave("FigA7_LibCon_Cont_CATE.jpeg", plot = combined_cate_plot, width=6 ,height = 7, dpi=600)

#######################################################################################################


#######################################################################################################
#### Study 3: Table A16: CATE Estimation, Continuous Operationalization of Lib/Con Ideology ####

exp_cate_models_short_label_c <- list("(EU)"=model4_c, "(UKR. Ref. Therm.)" =model1_c, "(Burden/Benefit)"=model2_c, "(Accept UKR Ref.)"=model3_c)

cate_coef_map <- c("treated"="Treatment", "lib_con"="Lib/Con", "treated:lib_con"="Treatment*Lib/Con")

modelsummary(exp_cate_models_short_label_c, output = "latex", coef_map = cate_coef_map, gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05))


#######################################################################################################

#######################################################################################################
#### Tables A17-A19 Mediation Models  ####

data$svk_control <- car::recode(data$svk_control , "1=0; 2=1; 3=2; 4=3; 5=4; 6=5; 7=6; 8=7; 9=8; 10=9; 11=10")
data$svk_eu_comp <- car::recode(data$svk_eu_comp , "1=0; 2=1; 3=2; 4=3; 5=4; 6=5; 7=6; 8=7; 9=8; 10=9; 11=10")
data$eu_control <- car::recode(data$eu_control , "1=0; 2=1; 3=2; 4=3; 5=4; 6=5; 7=6; 8=7; 9=8; 10=9; 11=10")


mediator_model_1 <- lm(svk_control ~ treated,
                     data = data)


outcome_model_a <- lm(ukr_therm ~ treated + svk_control,
                    data = data)

mediation_results_a <- mediate(mediator_model_1,   
                             outcome_model_a,      
                             treat = "treated", 
                             mediator = "svk_control",
                             boot = TRUE,        
                             sims = 1000)

summary(mediation_results_a)


outcome_model_b <- lm(burden_benefit ~ treated + svk_control,
                      data = data)

mediation_results_b <- mediate(mediator_model_1,   
                               outcome_model_b,      
                               treat = "treated", 
                               mediator = "svk_control",
                               boot = TRUE,        
                               sims = 1000)

summary(mediation_results_b)

outcome_model_c <- lm(ukr_accept ~ treated + svk_control,
                      data = data)

mediation_results_c <- mediate(mediator_model_1,   
                               outcome_model_c,      
                               treat = "treated", 
                               mediator = "svk_control",
                               boot = TRUE,        
                               sims = 1000)

summary(mediation_results_c)

outcome_model_d <- lm(eu_feel ~ treated + svk_control,
                      data = data)

mediation_results_d <- mediate(mediator_model_1,   
                               outcome_model_d,      
                               treat = "treated", 
                               mediator = "svk_control",
                               boot = TRUE,        
                               sims = 1000)

summary(mediation_results_d)

mediation_models1 <- list(
  "(EU)" = mediation_results_d,
  "(UKR. Ref. Therm.)" = mediation_results_a, 
  "(Burden/Benefit)" = mediation_results_b,
  "(Accept UKR Ref.)" = mediation_results_c
)


modelsummary(mediation_models1, 
             output = "latex", 
             stars = c('*' = .05),
             coef_omit = "acme_0|ade_0",
             coef_rename = c("acme_1" = "ACME Slovakia Control", 
                             "ade_1" = "ADE Treatment"))



mediator_model_2 <- lm(svk_eu_comp ~ treated,
                       data = data)


outcome_model_e <- lm(ukr_therm ~ treated + svk_eu_comp,
                      data = data)

mediation_results_e <- mediate(mediator_model_2,   
                               outcome_model_e,      
                               treat = "treated", 
                               mediator = "svk_eu_comp",
                               boot = TRUE,        
                               sims = 1000)

summary(mediation_results_e)


outcome_model_f <- lm(burden_benefit ~treated + svk_eu_comp,
                      data = data)


mediation_results_f <- mediate(mediator_model_2,   
                               outcome_model_f,      
                               treat = "treated", 
                               mediator = "svk_eu_comp",
                               boot = TRUE,        
                               sims = 1000)

summary(mediation_results_f)

outcome_model_g <- lm(ukr_accept ~ treated + svk_eu_comp,
                      data = data)

mediation_results_g <- mediate(mediator_model_2,   
                               outcome_model_g,      
                               treat = "treated", 
                               mediator = "svk_eu_comp",
                               boot = TRUE,        
                               sims = 1000)

summary(mediation_results_g)

outcome_model_h <- lm(eu_feel ~ treated + svk_eu_comp,
                      data = data)

mediation_results_h <- mediate(mediator_model_2,   
                               outcome_model_h,      
                               treat = "treated", 
                               mediator = "svk_eu_comp",
                               boot = TRUE,        
                               sims = 1000)

summary(mediation_results_h)

mediation_models2 <- list(
  "(EU)" = mediation_results_h,
  "(UKR. Ref. Therm.)" = mediation_results_e, 
  "(Burden/Benefit)" = mediation_results_f,
  "(Accept UKR Ref.)" = mediation_results_g
)

modelsummary(mediation_models2, 
             output = "latex", 
             stars = c('*' = .05),
             coef_omit = "acme_0|ade_0",
             coef_rename = c("acme_1" = "ACME EU Decision", 
                             "ade_1" = "ADE Treatment"))


mediator_model_3 <- lm(eu_control ~ treated,
                       data = data)


outcome_model_i <- lm(ukr_therm ~ treated + eu_control,
                      data = data)

mediation_results_i <- mediate(mediator_model_3,   
                               outcome_model_i,      
                               treat = "treated", 
                               mediator = "eu_control",
                               boot = TRUE,        
                               sims = 1000)

summary(mediation_results_i)


outcome_model_j <- lm(burden_benefit ~treated + eu_control,
                      data = data)


mediation_results_j <- mediate(mediator_model_3,   
                               outcome_model_j,      
                               treat = "treated", 
                               mediator = "eu_control",
                               boot = TRUE,        
                               sims = 1000)

summary(mediation_results_j)

outcome_model_k <- lm(ukr_accept ~ treated + eu_control,
                      data = data)

mediation_results_k <- mediate(mediator_model_3,   
                               outcome_model_k,      
                               treat = "treated", 
                               mediator = "eu_control",
                               boot = TRUE,        
                               sims = 1000)

summary(mediation_results_k)

outcome_model_l <- lm(eu_feel ~ treated + eu_control,
                      data = data)

mediation_results_l <- mediate(mediator_model_3,   
                               outcome_model_l,      
                               treat = "treated", 
                               mediator = "eu_control",
                               boot = TRUE,        
                               sims = 1000)

summary(mediation_results_l)

mediation_models3 <- list(
  "(EU)" = mediation_results_l,
  "(UKR. Ref. Therm.)" = mediation_results_i, 
  "(Burden/Benefit)" = mediation_results_j,
  "(Accept UKR Ref.)" = mediation_results_k
)

modelsummary(mediation_models3, 
             output = "latex", 
             stars = c('*' = .05),
             coef_omit = "acme_0|ade_0",
             coef_rename = c("acme_1" = "ACME EU Control", 
                             "ade_1" = "ADE Treatment"))
#######################################################################################################
